Sparse Representations for Structured Noise Filtering 

Bishnu P. Lamichhane* and Laura Rebollo-Neira* 
February 1, 2008 



Abstract 

The role of sparse representations in the context of structured noise filtering is dis- 
cussed. A strategy, especially conceived so as to address problems of an ill posed nature, 
is presented. The proposed approach revises and extends the Oblique Matching Pursuit 
technique. It is shown that, by working with an orthogonal projection of the signal to 
be filtered, it is possible to apply orthogonal matching pursuit like strategies in order to 
accomplish the required signal discrimination. 

1 Introduction 

The problem of structured noise filtering is introduced in [1] , where a number of relevant 
signal processing applications are discussed. It can be posed as follows: consider that 
a signal /, represented as an element of an inner product space Tt, is produced by the 
superposition of two components, f% and /2, each of which belongs to a different subspace 
of Tt. More precisely, / = f\ + fi with f\ 6 S\ C Tt and fa G S2 C Tt. Structured noise 
filtering (to be also termed signal discrimination or signal splitting) consists of singling 
out a particular component from the signal /. Provided that S\ and £2 are given, such 
that S\ n £2 = {0}, one component, say /1, can be extracted from / by an oblique 
projection onto S\ and along £2- On the contrary, the situation <Si PI £2 7^ {0} implies 
that the signal decomposition is not unique and the splitting can not be tackled in a 
straightforward manner by oblique projections. Moreover, even when theoretically the 
condition S\ n<S2 = {0} is satisfied, if the subspaces 1S1 and £2 are not well separated, the 
construction of the corresponding projector becomes ill posed. Consequently, the signal 
splitting can not be achieved by numerical calculations in finite precision arithmetics. 
Here we focus on such a situation. We assume that the given subspaces S\ and S2 are 
theoretically disjoint, but close enough to yield an ill posed problem. 

Our proposal for the numerical realization of the signal splitting is focussed on the 
search of a subspace of the given <Si, where a class of signals is considered to lie. It will be 
assumed throughout the paper that the class of signals to be considered is iT-sparse in a 
spanning set for S\. By this we mean that given a spanning set for S\ the corresponding 
linear superposition of a signal has at most K nonzero coefficients. The K- value should 
be less than or equal to the dimension of the subspace S r C S\ for which the construction 
of an oblique projection onto itself, and along S2, is well conditioned. This assumption is 
quite realistic, considering that in practice there is often a lack of complete knowledge on 
the actual subspace 5i and to be on the safe side one may overestimate it. 
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The main motivation of this paper is to highlight the essential role that sparse rep- 
resentations play in the problem of structured noise filtering. Such representations have 
been the subject of considerable work over the last ten years [2-14]. We will dedicate spe- 
cial attention to discuss and illustrate 'why' and 'how' sparse representations are relevant 
in the present context. 

A technique, termed Oblique Matching Pursuit (OBMP), has been recently advanced 
in relation to the above described problem [15]. Such a technique evolves by stepwise 
selection of the sought subspace. The selection criterion is based on the consistency 
principle [16,17]. In this communication we revise and extend the OBMP technique. We 
show that by working with a particular projection of the signal at hand, rather than with 
the signal itself, one can make use of previously proposed orthogonal matching pursuit 
like methodologies, so as to look for the signal subspace yielding the correct splitting. 

The paper is organized as follows: In Section [2] we introduce the mathematical setting 
for signal representation to be adopted here, together with a discussion on the construction 
of oblique projections. Section [3] highlights the importance of the search for sparse solu- 
tions in the construction of oblique projectors for structured noise filtering. The proposed 
strategy is discussed in Section 0J The conclusions are presented in Section [5j 

2 Mathematical Framework 

We consider a signal, /, to be an element of an inner product space ft. The square norm 
||/|| 2 is then induced by the inner product that we indicate as (/, /) and is defined in 
such a way that if a is a number, (af,f) = a*(f,f), with a* representing the complex 
conjugate of a. For the purpose of this contribution we assume that all the signals of 
interest belong to some finite dimensional subspace V of TL. Thus, there exists a finite set 
{vi € 7i}f£i spanning V. Consequently, for every signal in V there is a set of numbers 
{ c i\iL\ which allows us to express the signal as the linear superposition 

M 

f = ^CiVi, 

i=i 

which is also called atomic decomposition. 

Although a signal was defined as an element of an abstract inner product space, for 
processing purposes we need a numerical representation of such an object. The process 
of transforming a signal into a number is refereed to as measurement or sampling. The 
mathematical operation performing such a transformation is then a functional. Since 
considerations will be restricted to linear measurements, we represent them by linear 
Junctionals. Thus, making use of Riesz theorem [18] we can express a linear measurement 
as m = (w, f) for some w € H. Considering now M measurements mi, i = 1, . . . , M, each 
of which is obtained by a measurement vector Wi, we have a numerical representation of 
/ as given by 

m i = (w i J), i = l,...,M. (1) 

The question concerning the possibility of reconstructing / E V from measurements ob- 
tained with vectors in a different subspace has been addressed in [16, 17, 19, 20]. It is in 
principle obvious that every signal in V can be reconstructed from vectors {wi € TQf^ 
spanning a subspace WcM, provided that those vectors give rise to a representation of 
any projector onto V. The difference in using one projector or another appears when the 
projector acts on signals outside a subspace. We summarize next some features relevant 
to the construction of projectors. 
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2.1 Oblique projectors 

Every idempotent operator is a projector. Hence, an operator E is a projector if E 2 = E. 
The projection is along its null space and onto its range. When these subspaces are 
orthogonal E is called an orthogonal projector, which is the case if and only if E is 
self-adjoint. Otherwise it is called oblique projector. 

Given two closed subspaces, V € TL and W 1 - € 7i, such that S = V + W ± and 
V D W 1 - = {0}, the oblique projector operator onto V along W 1 - will be represented as 
E vw ±. Then E VV[! ± satisfies Ey W± = E vv ^± and, consequently, 

Ey W ±f = /, if feV 
E vw xf = 0, if /GW X . 

In the particular case for which W 1 " = V 1 - the operator E vv ± is an orthogonal projection 
onto V. For indicating an orthogonal projector onto a subspace, X say, we use the 
particular notation P%. 

Consider that {vi}f£ 1 is a spanning set for V and {u^fL^ a spanning set for W, which 
is the orthogonal complement of W 1 " in S, i.e., 5 = W © W -1 ", with © indicating the 
orthogonal sum. Thus the spanning sets of V and W satisfy = {Py\> v i\fLi- If the 

condition V D W -1 " = {0} is fulfilled, the operator E vw ± can be constructed as 

M 

E vw^ = ^2 v i( w h-), ( 2 ) 

i=\ 

where the operation (wi, •) indicates that E V yy± acts by performing inner products. The 
vectors {wi}^ in ([2]) are obtained from vectors {ui\fL 1 through the equation [17,21] 

M 

w i = J29i,j u i ( 3 ) 
j'=i 

with g| ■ the element of a matrix a pseudo inverse of the matrix G the elements 
of which are given by the inner products (ui,vj),i,j = 1, . . . , M. 

Remark 1. The pseudo inverse allows for the possibility of the spanning sets of V and 
W being redundant. However, the condition V n W -1 " = {0} implies that V and W should 
have the same dimension and therefore the rank of G equals the dimension of V and W. 

For later convenience, we introduce at this point an alternative representation of 
-Byyy±. To this end, denoting as ej,i = 1,...,M the standard orthonormal basis in 
C M , we define the operators V : C M -> V and W : C Af -> W as 

A/ M 
i=l i=l 

Thus the corresponding adjoint operators W* : W — > C M and V* : V — > C Ajf are 

M A/ 

i=l i=l 
It follows from the above definitions that Av^ = and #*Pyw = W*. 
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Considering that ip n € C M , n = 1, . . . , M, are the eigenvectors of matrix G = W*V = 
W*W, and assuming that there exist N nonzero eigenvalues A n , n = 1, . . . , N, on ordering 
these eigenvalues in descending order we can express the matrix elements of the Moore- 
Penrose pseudo inverse of G as: 

N 1 

^ = E^)ir<W' (4) 

n=l 

with ip n (i) the i-th component of ij) n . Moreover, the orthonormal vectors 

t ^ n /V 1 AT fZ\ 

U = , o- n = V^n, n=l,...,N (5) 

are singular vectors of W, which satisfies W*£ n = cr n i(i n , as it is immediate to verify. By 
defining now the vectors rj n , n = 1, . . . , N as 

^ n r ivr /vrt 

% = , n = l,...,N, (6) 

the projector E vv ^i_ in (J2j) is recast in the fashion 

N 



EvW ± =^2 r ln(€n,-)- (7) 

n=l 

Inversely, the representation ([2]) of -Eyyyi arises from ([7]) , since 

- 1 

i = l s ...,M. (8) 



AT 



n=l a ™ 

Proposition 1. T/ie vectors £ n S W, n = 1, . . . , N and r] n € V, n = 1,...,N given in 
([5]) and © are biorthogonal to each other and span W and V, respectively. 

Proof. Using ([5]) and ([6]) we have 

{£m,Vn) = (Wlp n ,Vlp m ) = (lp n ,W*Vlp m ) = 5 n ,m — = S n>m , (9) 



which proves the biorthogonality property. 

The proof that span{£ n }^ =1 = W stems from the fact that W = span{nj}^£ 1 = 
span{wj}^£ 1 , which allows us to express an arbitrary j £ Was the linear combination 
9 = Y,i=i a i w i- Then, using ©, we have g = I]^ = i a n £ n with d n = ^ £ i=1 aiVCWi 
which proves that W C spanj^}^. On the other hand for g G span-j^i}^^ we can 
write g = S n=1 dn£n and using © we have / = Yji=i dm, with d\ = i £ n=1 d n ip n (i). 
This proves that spanl^}^ C W and therefore span{£ n }^ =1 = W. The proof that 
span{T/ ra }^ =1 = V is equivalent to the previous one. □ 

Since E vw ±f = f for every signal / G V, regardless of the subspace W , one can 
consider a different subspace to construct measurement vectors and reconstruct a 
signal in V using any set of such vectors, as long as VnW^ = {0} with the orthogonal 
complement of W = span{wi}^£ 1 . On the other hand, if / ^ V, the measurement vectors 
can be chosen to be suitable for the particular processing task. For instance, if the goal 
is to produce an approximation fy of / in the subspace V, then in oder to minimize 
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the distance ||/ — /y|| we need /v = Pyf- Any other projection would yield a distance 
||/ — E V y V ±f\\ which satisfies [16] 



H/-iWll<ll/--EvW||< 



1 



Wf-PvfW 



cos(#) 



where 6 is the minimum angle between the subspaces V and W. The equality is attained 
for V = W, which correspond to the orthogonal projection. 

However, if the aim were to discriminate from a signal produced by different phenom- 
ena only the component in V, then, as discussed below, an oblique projection turns to be 
appropriate. 

Suppose that a signal / is the superposition of two signals f\ and fi with f\ € V 
and /2 in W -1- . The projection that will rescue f\ from / is E V y V ±f. A number of signal 
processing examples where an oblique projection is required are given in [1]. Provided 
that the subspaces hosting the signal components are well separated, the discrimination 
of components with different structure is successful. Unfortunately, this is not always the 
case and the construction of the necessary projector may generate an ill posed problem. 

3 The need for sparse representations in the present 
context 

This section is dedicated to illustrate, by recourse to a numerical example, the crucial 
role that the search for sparse solutions plays in the construction of oblique projectors for 
signal discrimination. Consider that the spaces V and W -1- , such that V R W 1 - = {0}, are 
given, and the spanning set for V is a basis of dimension M . For constructing the dual 
vectors Wi as in ([8]) we first construct Pyy± , to generate the vectors 



spanning W. 

Remark 2. Since VRW -1- = {0}, and the set {v is assumed to be linearly independent, 
the set {lij}^ is also linearly independent. 

Proof. Suppose that ^i u i = ® f° r some se t of numbers {bi}^ . Then (fTOj) implies 

g = Pyy±g for g = J2i=i ^>i v i- Since by definition g € V, and VnW 1 - = {0} by hypothesis, 
we conclude that g = 0. Hence, the fact that the set is linearly independent implies 

that b{ = 0, i = 1, . . . ,M, which establishes that the set {u^f^ is linearly independent. 



Remark 3. Conversely, the fact that nonzero vectors constructed as in (jlOp are linearly 
independent implies that V H = {0} [15]. 

In order to render the numerical calculation of the vectors Wi, i = 1, . . . , M spanning 
W as stable as possible, it is convenient to orthogonormalize vectors Uj, i = 1,...,M 
to obtain the vectors qi, i = 1, ...,M satisfying (qi,qj) = 5ij. With these vectors we 
construct the M x M matrix G having elements (qi, Vj). For the situation considered here 
this matrix has an inverse. Let us denote the element (i, j) of G -1 as g~j and construct the 
corresponding vectors Wi, i = 1, . . . , M as prescribed in (|3]) or ([8]). As will be illustrated 
by the numerical example below, in spite of the fact that 'theoretically' V Pi W -1 " = {0}, 
numerical errors, due to the existence of small singular values, may cause the failure to 
find the unique signal splitting that theoretically one should expect. 



m = Vi - P w ±Vi, i = 1, . . . , M 



(10) 



□ 
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Example 1. Let V be the cardinal cubic spline space with distance 0.065 between 
consecutive knots, on the interval [0, 10]. This is a subspace of dimension M = 163, which 
we span using a B-spline basis {Bi(x), x € [0, 10]}|£|. The background we wish to filter 
belongs to the subspace W 1 - spanned by the set of functions {yi{x) = (x + l) -0 - 05 *, x G 
[0, 10]}f£ 1 . Here the inner product is defined as (/,<?) = J ° f(x)*g(x) dx, and all the 
integrals are computed numerically. 

This example is very illustrative of how sensitive to numerical errors the computation 
of oblique projectors is. The subspace we are dealing with are disjoint: the last five 
singular values of the corresponding matrix G are: 

0.2305,0.2298,9.3211 x 10" 4 , 2.5829 x 10" 6 , 2.5673 x 10~ 7 , 

while the first is o\ = 1.5018. The smallest singular value cannot be considered a numerical 
representation of zero when the calculations are being carried out in double precision 
arithmetic. Hence, one can assert that the condition V Pi = {0} is fulfilled. However, 
due to the three small singular values the computation of the measurement vectors in 
the whole subspace W is inaccurate enough to cause the failure to correctly separate 
signals in V from their background. The left graph of Figure 1 is generated by a random 
superposition or K = 70 B-splines added to a background in the given W ± . The broken 
line in the right graph represents the oblique projection onto the given V along As 
can be seen, the projection does not produce the required signal, which is represented by 
the continuous dark line in the same graph. Now, since the spectrum of singular values has 
a clear jump (the last three singular values are far from the previous ones) it might seem 
that one could regularize the calculation by truncation of singular values. Nevertheless, 
such a methodology turns out to be not appropriate for the present problem, as it does 
not yield the correct separation. The light lines in the right graph of Figure 1 depict the 
three approximations obtained by neglecting one, two and three singular values. 

Propositions [2] below analyzes the effect that regularization by truncation of singular 
values produces in the resulting projection. 




Figure 1: Left graph: signal plus background. Right graph: the dark continuous line cor- 
responds to the signal to be discriminated from the one in the left graph. The broken line 
corresponds to the approximation resulting from the oblique projection. The light lines corre- 
spond to the approximations obtained by truncation of singular values (the one closest to the 
required signal correspond to the truncation of three singular values). 
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Proposition 2. Truncation of the expansion ([7]) to consider up to r terms, produces 
an oblique projector along VV r = W -1 + Wo + Vq, with Vq = span{r7j}^ r+1 and Wo = 
span{£i}£L r+1 , onto V r = span{r?j}[ =1 . 

Proof. The biorthogonality between {£}i=i an d {Vi}i=i established in Proposition [T] en- 
sures that = Yl r i= i m{£» •) is a projector, since = Ey^. 

As established in Proposition [Q V = spanj^j}^, and therefore every / 6 V can be 
decomposed as / = f r + f a with f r G span{r/j}[ =1 and f G span{r/j}^, r+1 . Moreover, 
Ey yfy f = f r , Ey yy f r = f r , and Eyy^f = 0, which proves that the projection is onto 
V r and Vo is included in the null space of Ey yy . Equivalently, for every g Q G Wo = 
span{^j}^ r+1 we have Ey yy g a = 0, because the set {Ci}iL\ is orthonormal. Thus, Wo is 
included in the null space of -E-y r yy r • □ 

3.1 Getting ready for a greedy search of the sparse solution 

We discuss here the properties that will be of assistance in the next section, where we 
will present our strategy for the search of the sparse representation achieving the desired 
signal discrimination. The goal is to avoid the computation of the measurement vectors 
in the whole subspace. Instead, we strive to find the subspace Vk C V, where the signal 
component one wants to discriminate from the noise is assumed to lie. We work under 
the hypothesis that the subspace W -1- is given and fixed. Furthermore, V n W ± = {0}, 
which implies that there exists a unique solution for the signal splitting. The problem 
we need to address arises from the fact that, if the subspaces V and W ± are not well 
separated, the numerical calculation of the measurement vectors is not accurate (due to 
the numerical operations being carried out in finite precision arithmetic). As a conse- 
quence, the representation of the corresponding projector fails to produce the correct 
signals separation. 

Assuming that we are able to accurately compute in finite precision arithmetic r 
measurement vectors, we could attempt to filter structured noise of a signal belonging 
to a subspace spanned by at most r of such vectors (i.e. the expansion of the signal in 
V should have at most r nonzero coefficients). However, even possessing this knowledge 
about the signal, the problem of finding the right subspace would be in general intractable: 
out of a set of cardinality M there exist (^) possible subsets of cardinality r. An adaptive 
strategy for the subspace selection, given a signal, is advanced in [15]. Before revising and 
extending that strategy we need to recall two relevant properties of oblique projectors. 

Property 1. The oblique projector E vyv ± satisfies Py\;Ey yv ± = Py^. 

Proof. It readily follows by applying Pyv on both sides of ([2]) or ([7]). Since (ui,Vj) = 
(ui ,Uj), one has 

M 

PyvEy W ± = ^UiiWi, ■) = P W . (11) 
t=l 

Moreover, since Pyv^i = £i> considering ([JJ we have 

N 

PwE vw x = ^e.fe, •) = Pw- (12) 
1=1 

□ 
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Property 2. Given a signal f inV + W -1- = W © W 1 - , the only vector g £ V satisfying 

Pwf = Pw9 (13) 

zs 5 = E vw ±f. 

Proof. If g = E V y^±f (]13p trivially follows from Property [H Let us assume now that 
there exists g G V such that (fl"3|) holds. Then Pyv(/ — 5) = 5 i-e., (/ — ff) £ W -1 -. Hence 
Ey W ±(f — g) = and, since g £ V, this implies that E V w_i_f = <?. □ 

Let us suppose that V& = span{tij}^ =1 is given and the spanning set is linearly inde- 
pendent. Assuming that Vk H W 1 - = {0} we guarantee that the set of vectors {ui}^ =1 , 
with m given in (I10p is also linearly independent. Therefore the dimension of Vk is equal 
to the dimension of Wk = span{uj}^ =1 = span{wf}f =1 . We use now a superscript k to 
indicate that the measurement vectors {wf }^ =1 span Wfc. Hence these vectors give rise to 
the oblique projection of a signal /, onto Vk and along W -1 , as given by: 

k k 
i=l i=l 

It is clear from ([14p that if the atoms in the atomic decomposition were to be changed (or 
some atoms were added to or deleted from the decomposition) the measurement vectors 
wf, and consequently the coefficients c\ in (|14p . would need to be modified. The recursive 
equations below provide an effective way of implementing the task. 

Forward /backward adapting of measurement vectors 

Starting with w\ = rr^rp, an d U\ as in (fT0|) . the measurement vectors wf +1 . 



1 . . . , k + 1 can be recursively constructed from wf, i = 1 . . . , k as follows [21]: 

w k + l = wf -w k k Xl(u k +i,wf), i = l,...,k (15) 

w kli = 71 TT2' Tfc+i =Wfe+l --fWfcWjfc+l! (16) 

ll7Jfe+l|| 

where Py^ k is the orthogonal projector onto Wk = span{«j}| =1 . We note that, since 
PyvWi = wf, (fT5|) can also be written as 

= w*-«/£g(t; fc+ i,u£), » = l,...,fc. (17) 

It follows from the above equations that when incorporating a linearly independent atom 
Vk+i in the atomic decomposition (|14p . the coefficients can be conveniently modified ac- 
cording to the recursive equations 



c 



ttl = «!»/>» ( 18 ) 
cf +1 = K fe+1 ,/) = ^-c^K fe ,^ +1 ), i = l,...,fc. (19) 

Conversely, considering that the atom, Vj say, is to be removed from the atomic decomposi- 
tion (fl"4"|) , and denoting the corresponding subspaces Vf.\j and Wf.\j , in order to span Wk\j 

k\ i 

the measurement vectors w i , i = 1, . . . , k are modified according to the equation [21] 

^ = vt- ^H\f\ i = l,...J-U + l,...,k. (20) 
\\ w j II 



Consequently, the coefficients in (|T4|) should be changed to 

^ = c?- ^. W V,^ ) , i i j I-./ • i /■•• (21) 

\\ w j\\ 
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4 Adaptive pursuit strategy for subspace selec- 
tion 



Given a signal /, we aim at finding the subspace Vk C V where the signal belongs. Let us 
stress once again that the problem arises from the impossibility of correctly computing the 
measurement vectors spanning the whole subspace W. Otherwise the right subspace is 
determined simply by the indices corresponding to the atoms having nonzero coefficients 
in the full atomic decomposition (I14j) . 

Since Vk+i + W = W^+i © W 1 - the forward selection criterion we propose is based 
on Property [21 which implies that if a given / satisfies Pyvf = P\v k f it also satisfies 
EyyyA-f = E VkW ±f. Thus, by fixing Py\) k , at iteration k + 1 we select the index £^+1 such 
that ||-Pyv/ — -P>V fe+ i/|| 2 is minimized. 



Proposition 3. Let us denote J to the set of indices i = 1, . . . , M. Given Wfc, the index 
4+1 corresponding the atom ug k+1 in the set {lijjiej f or which \\Pyyf — Py\) k+1 f\\ 2 *s 



minimal is to be determined as 

4+i = arg max V'^ 1 , 7 n + 0, (22) 

neJ\J k ||7n|| 

with 7 n given in (fTo]) , and the set of indices that have been previously chosen to deter- 
mine Wfc. 

Proof. It readily follows since Pw k+1 f = P\v k f + nj^lP^ anc ^ nence \\P\vf — Av fc+ i/|| 2 = 
\\Pwf\\ ~ \\PwJ\\ 2 - B ™se P w f and P W J are fixed, \\P w f - P Wk+1 f\\ 2 is 

minimized if np^j > 7n 7^ is maximal over all n G J \ . □ 

Remark 4. Since Pw k+1 = PwE Vk+lW ± we can write \\Pwf — Av fc+ i/|| 2 = \\Pw{f ~ 
E Vk+i y^± f)\\ 2 , and the condition of the previous proposition can be seen as the condition 
for minimizing the distance of E Vk+iW ±f to f , with respect to the weighted seminorm 
|| • \\p w induced by the weighted inner product (-,-)p w defined as {f,g)p w = {f,Pw9)- 

The OBMP selection criterion given in [15], which is based on the consistency principle 
[16, 17], selects the index 4+1 as the maximizer over n G J \ J*, of 

l<7n,/)| 1 1 I, / n 

-|| no - ' 7n f D. 

WlnV 

This condition was proposed in [15] so as to select the measurement vector w%t\ pro- 
ducing the maximum consistency error A = |(w^}, / — E Vk yy±f)\, with regard to a new 
measurement li^}- However, since the measurement vectors are not normalized to unity, 
it is sensible to consider the consistency error relative to the corresponding measurement 
vector norm Hu^t^H, and select the index so as to maximize over k + 1 G J\Jk the relative 
consistency error 

A= IK ^'f; + ^ W±/)l , HXlW^O (23) 
W w k+i\\ 

Property 3. The index 4+1 satisfying (|22p maximizes over k + 1 G J \J}~ the relative 
consistency error ()23p 
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Proof. Since for all vector w^t\ given in (fT6|) E v y V ±w k t\ = and = ||7fc+i|| 1 

we have 

X \(lk + lj)\ 

Hence, maximization of A over k + 1 G J \ </& is equivalent to ()22p . □ 

It is clear at this point that the forward selection of indices prescribed by proposition 
(|22l) is equivalent to selecting the indices by applying the Optimized Orthogonal Matching 
Pursuit (OOMP) [22] strategy on the projected signal Pyyf using the dictionary {itj}j g j. 

The hypothesis that the computation of more than r measurement vectors becomes an 
ill posed problem enforces the forward selection of indices to stop if iteration r is reached. 
Nevertheless, the fact that the signal is assumed to be -fT-sparse, with K < r, does not 
imply that before (or at) iteration r one will always find the correct subspace. The r-value 
just indicates that it is not possible to continue with the forward selection, because the 
computations would become inaccurate and unstable. Hence, if the right solution was 
not yet found, one needs to implement a strategy accounting for the fact that it is not 
feasible to compute more than r measurement vectors. An adequate procedure is achieved 
by means of the swapping-based refinement to the OOMP approach introduced in [23]. 
As discussed below, it consists of interchanging already selected atoms with nonselected 
ones. 

Consider that at iteration r the correct subspace has not appeared yet and the selected 
indices are labeled by the r indices £%,... ,£ r . In order to choose the label of the atom 
that minimizes the norm of the residual error as passing from approximation Pw r f to 
approximation Py\> r ^.f we should fix the index of the atom to be deleted, £j say, as the 
one for which the quantity 

(24) 



is minimized i = 1, . . . , r [23, 24]. 

The process of eliminating one atom from the atomic decomposition (I14p is called 
backward step while the process of adding one atom is called forward step. The forward 
selection criterion to choose the atom to replace the one eliminated in the previous step 
is accomplished by finding the index £i,i = 1, . . . , r for which the the functional 

e n = -^TpTp> with v n = u n - Pw rXj u n , \\u n \\ ^ (25) 



is maximized. In our framework, using (|20p . the projector Pyv r \j 1S computed as 



PW rX , = PWr 



v r\i rvr II r 1 12 

II^JII 

Since Py^ r and are available, the computation of the sequence v n in (|25p is a simple 
operation. 

As proposed in [23] the swapping of pairs of atoms is repeated until the swapping 
operation, if carried out, would not decrease the approximation error. The implementation 
details for an effective realization of this process are given in [23], and MATLAB codes 
are available at [25]. Since there is no guarantee that at the end of the swapping of pairs 
of atoms the correct subspace has been found, the process can continue by increasing the 
number of atoms the swapping involves. At the second stage, in line with [26] we propose 
the swapping to be realized by the combinations of two backward steps followed by two 
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forward steps, provided that the interchange of the two atoms improves the approximation 
error. If at the end of the second stage the right subspace has not yet been found, the 
number of atoms involved in the swapping is increased up to three and so on. Notice that 
if the number of atoms to be interchanged reaches the value r the whole process would 
repeat identically. This is avoided by initiating the new circle with a different initial atom. 
Although convergence cannot be guaranteed, the above specified hypothesis ensure that 
the algorithm will stop when the correct signal splitting has been found. At such a stage 
one has Pwf = Pw r f with W r spanned by the selected atoms ,£ r . If the order 

K of sparseness of the signal is less than r a number of r — K coefficients in the atomic 
decomposition 

r r 

1=1 8=1 

will have zero value. 

4.1 Examples 

Firstly we applied the proposed strategy to the numerical simulation of Example 1, which 
is a very simple test for our method and therefore in a run of 50 simulations we could 
produce the correct signals splitting at the stage involving forward selection only. 

Example 2. For this example we have used the same background as in Example 

1, but a dictionary of B-splines spanning the same space as the basis. The dictionary 
consists of functions of broader support than the basis functions for the same space, and 
the translation parameter is reduced (for more details on the construction of B-splines 
dictionaries see [27], MATLAB codes are available at [25]). In this case, the spectrum of 
singular values of matrix G decreases continuously, as shown in the top graphs of Figure 

2. Since it is difficult to decide on where to truncate the singular values, for the sake of 
comparison with the proposed technique we made a signal dependent truncation. This 
was achieved by setting the number Q of singular values to be considered so as to minimize 

WPwf-P^fW, 

where Wq indicates the subspace spanned by the first Q singular vectors of the operator 
W (c.f. ©) Neither in this case the regularization by truncation of singular values was 
successful. The result is depicted by the lighter line in the right bottom graph of Figure 
2. The dark line plots the sought signal. 

By means of the proposed strategy we were able to find the correct signal splitting 
in the 50 simulations we ran. Only in one of the cases a re-initialization took place. 
Increasing the number of atoms in the simulated atomic decomposition up to 80, we also 
found the correct splitting in all the cases. In this simulation the re-initialization stage 
occurred in the case of 5 signals, out of the 50 signal in the run. By increasing the 
number of atoms up to 90, re-initialization took place in the case of 8 signals. As could 
be expected, due to the singular values decay, by increasing the number of atoms up to 
100 in the atomic decomposition we started to observe instability in the calculations. 

Example 3. Here the signal space is spanned by M = 210 vectors in R L , with 
L = 1000, given as 
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Figure 2: Top graphs: spectrum of singular values of matrix G in Example 2. The right graph 
plots the line in the left graph, but in logarithmic scale. 

Bottom graphs: The one on the left shows the signal plus background. The dark line in the 
right one is the signal component to be discriminated from the signal in the left graph. The 
light line represents the approximation obtained by truncation of singular values. The proposed 
approach reproduces exactly the dark line. 



The space of the noise is spanned by the set 

{yi = ^350000-0.005^ j = 1) ... )L} 4£0 

We ran 50 simulations, keeping the noise fixed and considering a different realization of 
the signal, which was generated as a linear combination of 90 vectors taken randomly 
from the given spanning set. 

The spectrum of singular values is depicted in the top graphs of Figure 3. In this 
case the signal dependent criterion for truncation does achieve the correct signal splitting. 
The left bottom graph of Figure 3 shows one of the realizations of the signal plus noise in 
the simulation. The dark line in the right graph of Figure 3 plots the exact signal. The 
approximation obtained by truncation of singular values is plotted with a lighter line, 
which cannot be distinguished from the dark one in the scale of the figure. It is clear 
from this result that in this case the signal does not have a significant component in the 
subspace spanned by the neglected singular vectors. Similar results are obtained in all 
the other realizations in the simulation. The norm of the error in this case is 0.53, while 
the mean value of the error norm with respect to the 50 cases is 0.78. 
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By applying the proposed strategy for searching the sparse representation we found the 
exact solution in the 50 cases. Re-initialization was necessary in 5 of the 50 realizations. 




200 400 600 800 1000 200 400 600 800 1000 



Figure 3: Top graphs: spectrum of singular values of matrix G in Example 3. The right graph 
plots the line in the left graph, but in logarithmic scale. 

Bottom graphs: The one on the left shows the signal plus noise. The dark line in the right one 
is the signal to be separated from the one in the left graph. The approximation obtained by 
truncation of singular values is plotted by a lighter line, which cannot be distinguished from 
the other in the scale of the figure. The proposed approach reproduces the exact signal. 



5 Conclusion 

The role of sparse representations in the context of structured noise filtering has been dis- 
cussed. The discrimination of signal components is achieved by an oblique projection onto 
the right subspace. Considerations were restricted to those cases for which the signal sub- 
space and the noise subspace are theoretically complementary, but the construction of the 
dual basis for the whole signal subspace yields an ill posed problem, due to the calculations 
being carried out in finite precision arithmetic. It was shown by numerical simulations 
that, if the signal is sparse in a spanning set for the signal subspace, the required signal 
splitting may be achieved by means of adaptive techniques capable of searching for the 
required subspace while maintaining stability in the calculations. Although convergence 
of the proposed strategy for adaptive subspace search is not guaranteed, the method is 
capable to stop when the correct signal splitting is accomplished. 
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For the sake of comparison, an alternative regularization technique based on adaptive 
truncation of singular values was analyzed. The main disadvantage of such a technique 
lies in the fact that regularization is performed by a change of subspaces. Consequently, 
in general, the technique does not produce the required signal splitting. Moreover, even 
when a satisfactory splitting is attained (c.f. Example 3) the method does not provide an 
indication that this is so. Except for the very particular case in which the signal at hand 
has zero projection onto the subspace spanned by the disregarded singular vectors, the 
exact solution cannot be produced by this technique. On the contrary, the approach based 
on the search for the sparse representation is capable of producing the exact solution when 
the method stops. Thus, if the algorithm has converged, one can assert that the signal 
splitting is the required one. 
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